Crystallization Mechanism of Hard Sphere Glasses 
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In supercooled liquids, vitrification generally suppresses crystallization. Yet some glasses can still 
crystallize despite the arrest of diffusive motion. This ill-understood process may limit the stability 
of glasses, but its microscopic mechanism is not yet known. Here we present extensive computer 
simulations addressing the crystallization of monodisperse hard-sphere glasses at constant volume 
(as in a colloid experiment). Multiple crystalline patches appear without particles having to diffuse 
more than one diameter. As these patches grow, the mobility in neighbouring areas is enhanced, 
creating dynamic heterogeneity with positive feedback. The future crystallization pattern cannot 
be predicted from the coordinates alone: crystallization proceeds by a sequence of stochastic micro- 
nucleation events, correlated in space by emergent dynamic heterogeneity. 



A supercooled liquid is thermodynamically unstable 
towards crystallization, but may nonetheless survive for 
long periods. For moderate supercooling, crystalliza- 
tion is delayed by a large free-energy barrier separating 
the metastable fluid from the solid; thermal crossing of 
this barrier is followed by rapid deterministic growth of 
the resulting supercritical nucleus. This process is de- 
scribed fairly well by so-called 'classical nucleation the- 
ory' (CNT) For deep supercooling, the barrier is low 
and can be crossed by spatially diffuse collective motion 
sometimes compared to spinodal decomposition [2M3] • In 
spinodal decomposition, as in CNT, the late stages in- 
volve deterministic nonlinear amplification of order pa- 
rameter fluctuations [5[. At very deep supercooling, mo- 
tion can become so slow that the fluid's structure vitrifies, 
forming a glass, before crystallites can emerge. Though 
such structural arrest clearly inhibits crystal formation, 
it may not guarantee the glass to be long-lived. (Indeed, 
even glasses that appear permanent may unexpectedly 
crystallize in a process known as 'devitrification' jfj.) To 
reliably formulate stable glasses, a better mechanistic un- 
derstanding is needed of how glasses can crystallize de- 
spite their structural arrest. 

It is challenging experimentally to study particle-scale 
dynamics in atomic or molecular systems, due to the 
small size of their constituents. To overcome this prob- 
lem, colloidal systems (whose particles are visible with 
microscopy) have been widely used to explore the physics 
of glasses 0, Q and crystallization (9l4ll). Indeed, hard- 
sphere colloids (where density and supersaturation play 
the roles of inverse temperature and supercooling) are ar- 
guably the simplest system displaying fluid-solid coexis- 
tence and glass formation. T hey have therefore inspired 
many computer simulations 12 Most of these ad- 



dress nucleation from the ergodic fluid; here CNT broadly 
holds although open questions remain Few sim- 

ulation studies have addressed crystallization of glasses 
— although this has been reported experimentally for col- 
loids @, [3] (as well as molecules p) — and no micro- 
scopic description of its mechanism yet exists. 

Here we explore the mechanism of glass crystalliza- 



tion (at fixed volume, as in colloidal experiments) at the 
level of individual crystallites and particles. We ana- 
lyze the emergence of crystalline patches (of which none 
are initially present), finding that these appear and grow 
without needing even a subset of particles to move much 
beyond a diameter (a stronger result than in 0). The 
growth of crystalline patches enhances the mobility of 
the surrounding particles, creating an emergent dynamic 
heterogeneity (DH), not present initially, which gives pos- 
itive feedback for further crystallization. Notably, even 
in the late stages, domain growth is not controlled by 
barrier- free (quasi-deterministic [5]) evolution of any con- 
figurational order parameter. Instead, the future pat- 
tern of crystallinity depends on thermally random parti- 
cle velocities; reassigning these, even midway through the 
domain growth, leads to a different final pattern. This 
sensitivity resembles the stochastic, activated dynamics 
of early-stage CNT. Thus our crystallization mechanism 
for glasses can be viewed as a chaotic sequence of random 
micro- nucleation events, correlated in space by emergent 
DH. (To determine whether this process is related to the 
putative 'spinodal crystallization' of non-vitrified fluids 
at deep undercooling [l-Q, Hq| requires detailed clarifica- 
tion of the latter mechanism.) 

In the results presented here we use as units the par- 
ticle mass m, diameter a, and the thermal timescale 
t = a(m/k B T)^ 2 (so that m = a = k B T = 1). Our 
event-driven molecular dynamics (MD) study addresses 
N = 86400 monodisperse hard spheres in an NVT en- 
semble, with periodic boundary conditions in a cubic 
box of volume V. (For more details of the methods see 
[HI, Hll, Hll ) • To analyze the results, each particle is first 
assigned a vectorial bond order parameter d§ [23|. Pairs 
of neighbours are then identified using the criterion of 
van Meel and Frenkel [24J, and a pair is deemed 'con- 
nected' if the scalar product of their d§ vectors exceeds 
0.7. (The prescription of [24| avoids a dependence on an 
arbitrary cutoff distance used in earlier work [22]). Each 
particle is then labeled as 'crystalline' if it is connected 
with at least 6 neighbours, and as 'amorphous' other- 
wise; the fraction of crystalline particles is then denoted 
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X(t). We further quantify the local degree of crystalline 
order around any particle through a scalar bond-order 
parameter Qq as defined in Ref. [25j. 

In order to generate a glassy initial configuration with- 
out pre-existing crystallization nuclei, we compressed a 
small system (400 particles) to a high packing fraction, 
(j> w 0.64, and then replicate it periodically in space af- 
ter checking that the fraction of solid-like particles is 
less than 0.005. The configuration thus generated was 
isotropically expanded to the desired density before start- 
ing the MD run. In this way make sure there are no 
non-amorphous regions (Qq > 0.25 EBI ) m the initial 
configuration. Note that after a short transient, no dis- 
cernible trace of the initial periodicity could be detected 
(Fig. |4]). Periodic boundaries do not affect the appear- 
ance or the growth of crystallites provided that the sys- 
tem is large enough to accommodate a good number of 
clusters (Fig. 1(b)). 




FIG. 1: (a) fraction of crystalline particles (X) versus time 
for single runs at several packing fractions (see legend). The 
difference between curves 0.61 a and b is discussed in Fig. 4. 
(b) growth of the largest crystalline cluster for cj) = 0.54 (top) 
and 4> = 0.61 (bottom). 

Fig. dja) shows the fraction of crystalline particles, 
X(t), for various volume fractions 0.54 < <p < 0.63. 
At <p = 0.54 the fluid is metastable for a short period 
(30 < t < 60) until a post-critical cluster nucleates and 
grows rapidly after t ~ 100. For <j> > 0.55, X(t) starts 
growing almost immediately {t ~ £ = 1): hence no sig- 
nificant nucleation barrier is present. Within the ergodic 



fluid at 0.55 < <j> < 0.58, X(t) barely depends on (2fi|, 
while in contrast, for glasses (0 > <j> g ~ 0.58), X(t) grows 
more slowly initially, but increases its slope sharply once 
sufficient particles have crystallized to increase the free 
volume available to the remainder. Fig. Wfo) contrasts 
the domain growth pattern of crystallites in the CNT 
regime ((f) = 0.54) with that in glass crystallization. In- 
stead of a growing compact nucleus, a branched crys- 
talline network develops. This is shown for <fi = 0.61, the 
highest (f) for which we observed crystallization within 
our time window (Fig. Ufa)). From now on we focus on 
this case. 

In [l9|| we established that monodisperse hard sphere 
glasses can crystallize, even without pre-existing nuclei, 
while the root-mean-square (rms) particle displacement 
remains less than a particle diameter. Thus crystalliza- 
tion occurs without macroscopic diffusion — which is 
why it remains possible within the glass. But that anal- 
ysis did not rule out a form of DH involving larger dif- 
fusive displacements of the minority of particles actually 
involved in early-stage crystallization. This is however 
ruled out by Fig. [2j which shows separate displacement 
distributions for particles that become crystalline during 
a given time window and particles that remain amor- 
phous: most crystallizing particles moved far less than a 
diameter since the quench. Fig. [2] (inset) shows, for each 
population, the probability distribution F with which a 
fraction f c of a particle's neighbours were also neighbours 
at the start of the run. Even for crystallizing particles, 
(fc) > 0-7. Thus crystallites indeed form without large 
diffusive rearrangements (cage-breaking), requiring only 
local shuffling of particles from amorphous into ordered 
patterns. Visual observation of particle trajectories con- 
firms this [I?}. Nonetheless, the mean-square displace- 
ments of crystallizing particles are measurably larger, 
and the mean f c smaller, than those of their amorphous 
counterparts (Fig. [2]). This suggests a more subtle form 
of DH may contribute to the crystallization mechanism. 

To elucidate its role, we present in Fig.[3]maps showing 
the crystalline particles and the 5% most mobile parti- 
cles. While there is little DH initially, the most mobile 
particles become highly clustered once crystallites start 
to form, are preferentially found next to these crystallites, 
and have a higher tendency to become crystalline there- 
after. The dynamics not only involves slow growth of 
crystallites but also the preferential creation of new ones 
next to the old, in regions of enhanced mobility (Fig. [3] 
B — > C). The crystallites are denser than their surround- 
ings, as would be expected if pressures are fairly uniform 
(e.g., at X = 0.08 crystallites have <p c± 0.65 well above 
the global mean </> = 0.61). Thus our emergent DH is 
probably the result of an increase in free volume locally. 
Recently Ref. [Hj observed that particles at the inter- 
face between ordered and amorphous regions have lower 
density than average, although no feedback between DH 
and ordering emerged from that work. 
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FIG. 2: Probability distributions of displacements between 
the start of the run and time to.i defined by X(to.i) =0.1. 
Two distributions are shown: in black are particles that were 
amorphous at t = to.os and which remain amorphous at to.i. 
In red are particles that were amorphous t = to.os but are 
crystalline at to.i. (Note that once a particle crystallizes it 
barely moves, which is why we examine only the mobility of 
particles that crystallized near the end of the selected time 
window at to.i.) As can be seen in F(Ar), most of the parti- 
cles that crystallize do not move beyond one particle diameter 
during crystallization. Particles in the tail of the distribution 
do not have a crucial role as many crystalline regions appear 
without any of these being involved. Inset: probability distri- 
bution of the fraction f c of a particle's neighbours at time to.i 
that were also neighbours at t = ('conserved neighbours'). 
Colour code is the same as in main figure. 



Two runs starting from the same initial particle coordi- 
nates but with different (Maxwellian) initial velocities are 
shown in Fig.SJ The domain patterns (measured by Qq) 
evolve quite differently. This confirms the absence of pre- 
existing nuclei in our initial configuration: crystallinity 
arises at random locations. Indeed our initial coordi- 
nates, but not velocities, have a higher periodicity than 
the simulation box, which the emerging crystallinity pat- 
tern clearly ignores. Moreover, even if the particle veloc- 
ities are re-assigned midway through the domain growth, 
the subsequent evolution is markedly altered - including 
that of crystallites that already exist (Fig. [4]). There are 
an infinite number of possible chaotic sequences of parti- 
cle collisions in a given region; it appears that only some 
of these allow the crystalline order to increase. Thus the 
fate of a given region can be switched from amorphous 
to crystalline by changing either the initial or the mid- 
run velocities. The stochastic nature of the process by 
which particles in the glass incorporate to growing crys- 
tallites explains the difference between curves 0.61 a and 
b in Fig. la. As stated previously, this chaotic dynamics 
rules out any late-stage spinodal-like description based 
on continuous free energy reduction by evolution of a 
configurational order parameter. 

Examples of such order parameters are not only the 






FIG. 3: Slab in the xy plane showing the 5% most mobile 
particles (in red) and the crystalline particles (in light blue) 
at time t = (A), t = 320 (B), t = 640 (C) and t 1280 
(D). Mobile particles are ranked by the distance they move 
between the time of the frame at which they are shown and 
the subsequent frame (accordingly, frame D does not show 
mobile particles). Mobile particles are spatially correlated 
with crystalline ones and have a higher tendency to become 
crystalline than "average" amorphous particles. 



density but also orient ational parameters, including de 
and Qq, as used to identify regions of 'medium-range 
crystalline order' (MRCO). The role of such configura- 
tional order parameters in supercooled liquid dynamics 
is much debated [3]. For ergodic supercooled fluids, a 
two-step mechanism by which crystals nucleate within 
MRCO regions is increasingly accepted TH-l?)- For crys- 
tallization from the glass, Fig. |4j high-Qe regions seen at 
late times likewise correspond to intermediate- Q 6 ones 
at earlier stages (though in our simulations, none are 
present initially). However this conversion is quite grad- 
ual and does not involve an identifiable nucleation step. 
In [17j it was further argued that diffusionless crystal- 
lization might arise via "positional ordering in a region 
already having hep-like bond orientational order" . Such 
regions need not be pre-existing, however: indeed all are 
eliminated from our initial conditions. A concept recently 
introduced to address DH in glasses is that of 'propensity' 



28|, |29| . The basic idea of propensity is that configura- 
tional order does not directly determine the motion of a 
neighbourhood, but influences its probability of under- 
going motion. In this language an order parameter such 
as Qq could create a propensity for crystallization. 

Above we studied only 'fresh' glasses, immediately 
post-quench, whose amorphous structure undergoes ag- 
ing until crystallization intervenes. Using methods devel- 
oped in [26| , in future work we hope to find whether de- 





X= 0.002 



FIG. 4: Slab of the system for various values of the evolv- 
ing crystallinity X. Particles are coloured according to the 
degree of crystalline order in their neighbourhood (blue, < 
Q Q < 0.15; green 0.15 < Q 6 < 0.25; yellow 0.25 < Q 6 < 0.35; 
orange 0.35 < Q 6 < 0.45; red 0.45 < Q 6 < 0.55). The initial 
state (left) has a periodic density pattern that is quickly for- 
gotten. Its evolution follows two different trajectories (white 
arrows) from identical initial particle coordinates, but with 
different initial velocity choices (drawn at random from the 
thermal distribution). These trajectories correspond to curves 
labelled as <j) = 0.61a and b in Fig. 1. The black arrow shows 
a run where velocities are randomized after X reaches 0.05. 
The subsequent evolution is again altered, even though sig- 
nificant crystallinity was already present. At no stage do we 
find the future evolution to depend reproducibly on coordi- 
nates alone, although crystallites are more likely to form in 
regions of high Qq than elsewhere. 



vitrification of aged glasses follows the same microscopic 
mechanism as described here for fresh ones. 

In summary, we have elucidated the mechanism of 
crystallization of monodisperse hard-sphere glasses. As 
previously surmised [19] amorphous regions gradually 
transform into crystallites by 'shuffling' motion. Crys- 
talline regions grow in an autocatalytic way: shuffles 
leading to crystallization create mobile regions within 
which further crystallization is facilitated, giving an 
evolving dynamic heterogeneity that was not present ini- 
tially. The local emergence of crystallinity is chaotic: 
its future course depends on the chance values of parti- 
cle momenta, which determine which of many possible 
collisional trajectories actually arises. Structural order 
parameters do not fully determine where crystallites will 
form or how they will grow. The mechanism could be ex- 
perimentally validated via particle tracking in colloidal 
suspensions provided that the polidispersity allows for 
crystallization without fractionation (up to 3-4%). 
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